Thiele's interpolation formula

In mathematics, Thiele's interpolation formula is an interpolation formula for a function f(•) of a single variable, due to the Danish mathematician Thorvald Nicolai Thiele. It is expressed as a continued fraction, where ρ represents the reciprocal difference:

f(x) = f(x_1) %2B \cfrac{x-x_1}{\rho(x_1,x_2) %2B \cfrac{x-x_2}{\rho_2(x_1,x_2,x_3) - f(x_1) %2B \cfrac{x-x_3}{\rho_3(x_1,x_2,x_3,x_4) - \rho(x_1,x_2) %2B \cdots}}}

Here is a version of Thiele interpolation in Algol 68:

¢ The MODE of lx and ly here should really be a UNION of "something REAL", "something COMPLex",
and "something SYMBOLIC" ... ¢
PROC thiele:=([]REAL lx,ly, REAL x) REAL:
BEGIN
  []REAL xx=lx[@1],yy=ly[@1];
  INT n=UPB xx;
  IF UPB yy=n THEN
¢ Assuming that the values of xx are distinct ... ¢
    [0:n-1,1:n]REAL p;
    p[0,]:=yy[];
    FOR i TO n-1 DO p[1,i]:=(xx[i]-xx[1+i])/(p[0,i]-p[0,1+i]) OD;
    FOR i FROM 2 TO n-1 DO
      FOR j TO n-i DO
        p[i,j]:=(xx[j]-xx[j+i])/(p[i-1,j]-p[i-1,j+1])+p[i-2,j+1]
      OD
    OD;
    REAL a:=0;
    FOR i FROM n-1 BY -1 TO 2 DO a:=(x-xx[i])/(p[i,1]-p[i-2,1]+a) OD;
    yy[1]+(x-xx[1])/(p[1,1]+a)
  ELSE
    error; ~ ¢ Unequal length arrays supplied ¢
  FI
END;

Note that, although it works in most cases, it is sensitive to input values, especially those due to equally spaced abscissæ. (Essentially, in such a case, it may reduce to the ratio of two polynomials, which may have factors in common - thus yielding a 0/0 form.)

References